/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */

package uk.ac.liv;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileReader;
import java.io.IOException;
import java.util.*;

/**
 *
 * @author jonesar
 * Code assumptions/variables that can be altered:
 *
 * 1. Peptide can crosslink with itself - see peptideCanSelfCrosslink
 * 2. int minimumPeptideLength = 5;
 */
public class ProcessFasta {

    //Ultimately this could be parameterized for different enzymes
    private String enzymeRegex = "(?<=[KR])(?!P)";                  //defaults to trypsin
    //Some variables that affect the running of the algorithm
    private Boolean noCrossLinksOnCTerminus = true;
    private Boolean noCrossLinksOnModifiedRes = true;
    private Boolean peptideCanSelfCrosslink = true;
    private int minimumPeptideLength = 5;
    private String accessionRegex = " ";                            //assumption that the accession in the fasta file is followed by a space
    private int missedCleavages = 2;
    private HashMap<String, String> accToSeq = new HashMap();
    private HashMap<String, String> accToDefline = new HashMap();
    private HashMap<String, String[]> accToPeptides = new HashMap();
    private ArrayList<String> allProteinIDs = new ArrayList();
    private HashMap<String, Double> aaMap = new HashMap();
    private Double crossLinkerMass;       //These need to be parameterized and sent by the calling code
    private String crossLinkedResidue;       //These need to be parameterized and sent by the calling code 
    private static Double HMASS = 1.007825;
    private static Double OMASS = 15.994915;
    private static Double CTERM = 17.002740;
    private static Double NTERM = 1.007825;
    private HashMap<String, Double> varMods;    //hash map of var mods to be applied  ;key = "one letter code amino acids plus [ or ] for termini" value = double mass shift
    private HashMap<String, Double> fixedMods;  //hash map of fixed mods to be applied ;key = "one letter code amino acids plus [ or ] for termini" value = double mass shift
    private HashMap<String, Integer> varModStringToUniqueID = new HashMap();    //hashmap storing a unique ID for each var mod
    private HashMap<Integer, String> varModUniqueIDToString = new HashMap();    //reverse mapping

    /*
     * Constructor
     */
    public ProcessFasta(String fastaFile, Double crosslinkMass,
                        String crossLinkResidue, int missCleave,
                        HashMap varModifications, HashMap fixedModifications,
                        Boolean noCrossLinksOnCTerminalResidue,
                        Boolean noCrossLinksOnModifiedResidue) {


        crossLinkerMass = crosslinkMass;
        crossLinkedResidue = crossLinkResidue;
        missedCleavages = missCleave;
        varMods = varModifications;

        noCrossLinksOnCTerminus = noCrossLinksOnCTerminalResidue;
        noCrossLinksOnModifiedRes = noCrossLinksOnModifiedResidue;

        //Populate the hashmap giving unique IDs for all variable mods
        Integer id = 1;
        for (String varMod : varMods.keySet()) {
            varModStringToUniqueID.put(varMod, id);
            varModUniqueIDToString.put(id, varMod);
            id++;
        }

        fixedMods = fixedModifications;
        init(fastaFile);
    }

    private void init(String fastaFile) {
        fillAAMap();

        //Handle fixed mods
        for (String fixedMod : fixedMods.keySet()) {
            addFixedModsToAAMap(fixedMod, fixedMods.get(fixedMod));
        }

        readFasta(fastaFile);

    }

    private void readFasta(String inputFasta) {

        try {

            BufferedReader in = new BufferedReader(new FileReader(inputFasta));
            String line;

            String currSequence = "";
            String currProtAcc = null;
            String currDefline = null;

            int recordCounter = 0;
            while ((line = in.readLine()) != null) {
                line = line.replaceAll("\n", "");
                line = line.replaceAll("\r", "");

                if (line.contains(">")) {
                    //Insert previous into hash and reset
                    if (recordCounter != 0) {
                        currSequence = currSequence.replaceAll(" ", "");
                        allProteinIDs.add(currProtAcc);

                        String[] peptides = currSequence.split(enzymeRegex);
                        addPeptidesToMap(currProtAcc, peptides);

                        accToSeq.put(currProtAcc, currSequence);
                        //System.out.println("Inserting:" + currProtAcc + "_" + currSequence);
                        accToDefline.put(currProtAcc, currDefline);
                        //System.out.println("Inserting2:" + currProtAcc + "_" + currDefline);

                        currSequence = "";
                    }

                    line = line.replaceAll(">", "");

                    int splitPos = line.indexOf(accessionRegex);
                    if (splitPos != -1) {
                        currProtAcc = line.substring(0, splitPos);
                        currDefline = line.substring(splitPos + 1);
                    }
                    else {
                        System.out.println("Regular expression not found for split: " + line);
                        System.exit(0);
                    }

                    recordCounter++;
                }
                else {
                    currSequence += line;
                }
            }
            //handle last
            currSequence = currSequence.replaceAll(" ", "");
            accToSeq.put(currProtAcc, currSequence);
            accToDefline.put(currProtAcc, currDefline);

            String[] peptides = currSequence.split(enzymeRegex);
            addPeptidesToMap(currProtAcc, peptides);
            //System.out.println("Put last:" + currProtAcc + " seq: " + currSequence);
            allProteinIDs.add(currProtAcc);

            in.close();
        }
        catch (IOException e) {
            System.err.println("Error: " + e.getMessage());
        }
    }

    private void addPeptidesToMap(String protAcc, String[] peptides) {

        ArrayList<String> tempPeptides = new ArrayList();

        if (missedCleavages == 0) {
            accToPeptides.put(protAcc, peptides);
        }
        else {

            for (int i = 0; i < peptides.length; i++) {
                String outerPep = peptides[i];
                tempPeptides.add(outerPep);

                //TODO - test missed cleaves code!
                String tempPeptide = "";
                for (int j = i + 1; j <= i + missedCleavages; j++) {

                    if (j < peptides.length) {
                        tempPeptide += peptides[j];
                        String mc = outerPep + tempPeptide;
                        tempPeptides.add(mc);
                        //System.out.println("Adding mc:" + mc );
                    }
                }
            }

            String[] newPeptides = new String[tempPeptides.size()];
            for (int i = 0; i < tempPeptides.size(); i++) {
                newPeptides[i] = tempPeptides.get(i);
            }

            accToPeptides.put(protAcc, newPeptides);

        }

    }

    private void fillAAMap() {

        aaMap.put("A", 71.037114);
        aaMap.put("R", 156.101111);
        aaMap.put("N", 114.042927);
        aaMap.put("D", 115.026943);
        aaMap.put("C", 103.009185);
        aaMap.put("E", 129.042593);
        aaMap.put("Q", 128.058578);
        aaMap.put("G", 57.021464);
        aaMap.put("H", 137.058912);
        aaMap.put("I", 113.084064);
        aaMap.put("L", 113.084064);
        aaMap.put("K", 128.094963);
        aaMap.put("M", 131.040485);
        aaMap.put("F", 147.068414);
        aaMap.put("P", 97.052764);
        aaMap.put("S", 87.032028);
        aaMap.put("T", 101.047679);
        aaMap.put("U", 150.95363);
        aaMap.put("W", 186.079313);
        aaMap.put("Y", 163.06332);
        aaMap.put("V", 99.068414);
        aaMap.put("[", 0.0);             //N-terminus can be inserted for mods
        aaMap.put("]", 0.0);             //C-terminus (see getMR method)
        /*
         * A=71.037114
         * B=114.534940
         * C=160.030649
         * D=115.026943
         * E=129.042593
         * F=147.068414
         * G=57.021464
         * H=137.058912
         * I=113.084064
         * J=0.000000
         * K=128.094963
         * L=113.084064
         * M=131.040485
         * N=114.042927
         * O=0.000000
         * P=97.052764
         * Q=128.058578
         * R=156.101111
         * S=87.032028
         * T=101.047679
         * U=150.953630
         * V=99.068414
         * W=186.079313
         * X=111.000000
         * Y=163.063329
         * Z=128.550590
         * Hydrogen=1.007825
         * Carbon=12.000000
         * Nitrogen=14.003074
         * Oxygen=15.994915
         * Electron=0.000549
         * C_term=17.002740
         * N_term=1.007825
         */

    }

    /*
     * Helper method to change the mass of amino acid in the main hash i.e. due to a fixed modification
     */
    private void addFixedModsToAAMap(String aa, double modMass) {

        double mass = aaMap.get(aa);
        double newMass = mass + modMass;
        aaMap.put(aa, newMass);
        System.out.println("Changed mass of " + aa + " to " + newMass);
    }

    /*
     * Method for calculating the mass of a peptide plus a modstring e.g. PEPTIMER_0000001000 for Oxidation on M
     * where 1 is unique id for this mod. modstring is peptide length + 2 for N and C-terminus
     *
     */
    private double getPepMr(String peptidePlusModString) {
        double mr = 0.0;

        peptidePlusModString = peptidePlusModString.toUpperCase();
        peptidePlusModString = peptidePlusModString.trim();
        //System.out.println("Peptide:" + peptidePlusModString);

        String[] temp = peptidePlusModString.split("_");
        String peptide = "[" + temp[0] + "]";
        String modString = temp[1];

        String[] aas = peptide.split("");
        for (String aa : aas) {
            if (aa.length() == 1) {
                if (aaMap.get(aa) != null) {
                    mr += aaMap.get(aa);
                }
                else {
                    System.out.print("Non-recognized character in peptide:" + peptide);
                }
            }
        }
        String[] mods;
        //Code still functional if no modstring is provided
        if (modString != null) {
            mods = modString.split("");
            //Now loop through the var mod string
            for (String mod : mods) {
                if (!mod.equals("0") && !mod.equals("") && !mod.equals("X")) {
                    String varMod = varModUniqueIDToString.get(Integer.parseInt(mod));
                    //System.out.print("mod" + mod + " varMmod:" + varMod);
                    Double varModMass = varMods.get(varMod);
                    //System.out.println(" modmass:" + varModMass);
                    mr += varModMass;
                }
            }
        }


        mr += 2 * HMASS + OMASS;
        //System.out.println( " mr: " + mr); 

        return mr;

    }

    /*
     * Returns a hashmap of all peptide sequences (including potentially crosslinked), along with their mass
     * HashMap String = peptide sequence (plus special characters for crosslinker)
     * HashMap Double = corresponding mass
     *
     */
    public HashMap<String, Double> getAllPeptideMasses() {

        HashMap<String, Double> pepMasses = new HashMap();

        /*
         * All proteinIDs is an arraylist of all proteins seen (tested with only one protein so far
         * Use this to retrieve each array of peptides in turn
         *
         * Code as written does not search for crosslinked peptides in different proteins, only for peptides within one protein
         *
         */
        for (String protAcc : allProteinIDs) {

            //System.out.println("Processing: " + protAcc);           
            String[] peptides = accToPeptides.get(protAcc);


            for (int i = 0; i < peptides.length; i++) {
                String pep = peptides[i];
                if (pep.length() >= minimumPeptideLength) {

                    String modString = makeUnmodifiedModString(pep);  //Add the empty mod string 
                    String pepPlusModString = pep + "_" + modString;

                    Double pepMass = getPepMr(pepPlusModString);     //GetPepMR method is okay to receive the modstring after the peptide
                    pepMasses.put(pepPlusModString, pepMass);     //Add unmodified peptide to the hash
                    //System.out.println("Unmodified pep" + pep );
                    //System.out.println("\tAdd:" + pepPlusModString + " mr: " + pepMass);

                    //Handle var mods               
                    for (String varMod : varMods.keySet()) {

                        Double varModMass = varMods.get(varMod);
                        Set<Integer> modPositions = new HashSet();

                        int pos = pep.indexOf(varMod);
                        while (pos != -1) {
                            modPositions.add(pos);
                            pos = pep.indexOf(varMod, pos + 1);
                        }

                        //Generate x * y combinations of var mods via Powerset method
                        Set<Set<Integer>> modSets = SystemUtils.powerset(modPositions);
                        for (Set<Integer> modSet : modSets) {
                            if (modSet.size() > 0) { //not interested in empty set
                                //System.out.println("Processing set: " + modSet);
                                String newModString = insertModIntoModString(modString, varMod, modSet);
                                String modifiedPepPlusModString = pep + "_" + newModString;
                                pepMasses.put(modifiedPepPlusModString, getPepMr(modifiedPepPlusModString));
                                //System.out.println("\tAdd:" + modifiedPepPlusModString + " mr: " + getPepMr(modifiedPepPlusModString));
                            }
                        }
                    }
                }
            }
        }


        //Now this code only adds potentially crosslinked peptides to the hashmap if they both contain the crosslinked residue (unmodifed)
        Object[] peptides = pepMasses.keySet().toArray();
        for (int i = 0; i < peptides.length; i++) {

            String outerPepPlusModString = (String) peptides[i];
            String[] temp = outerPepPlusModString.split("_");
            String outerPep = temp[0];
            String outerPepModString = temp[1];
            Double outerPepMass = pepMasses.get(outerPepPlusModString);
            ArrayList<Integer> xlOuterPositions = new ArrayList();
            int pos = outerPep.indexOf(crossLinkedResidue);

            while (pos != -1) {

                //System.out.println(outerPep + " " + outerPepPlusModString + " " + noCrossLinksOnCTerminus + " " + pos + " " + (outerPep.length()-1) + " " + noCrossLinksOnModifiedRes + " " + outerPepModString.subSequence(pos+1,pos+2));
                //This code serves two purposes:
                //1) precludes crosslinks on C-terminal amino acid
                //2) then allows the check for the position not being modified
                if ((noCrossLinksOnCTerminus && pos != outerPep.length() - 1) && (noCrossLinksOnModifiedRes && outerPepModString.subSequence(pos + 1, pos + 2).equals("0"))) {
                    //System.out.println("\tAdd");
                    xlOuterPositions.add(pos + 1);    //Add position plus one, since the modstring has the n-terminus character at the beginning
                }
                else {
                    //System.out.println("\tDon't Add");                    
                }
                pos = outerPep.indexOf(crossLinkedResidue, pos + 1);
            }

            //Generate x * y combinations of potential crosslinked sites
            for (int x = 0; x < xlOuterPositions.size(); x++) {
                int xlPos = xlOuterPositions.get(x);
                String newOuterPepModString = outerPepModString.subSequence(0, xlPos) + "X" + outerPepModString.subSequence(xlPos + 1, outerPepModString.length());
                outerPepPlusModString = outerPep + "_" + newOuterPepModString;
                //System.out.println("new outer X:" + outerPepPlusModString);
                //String newOuterPep = outerPep.subSequence(0, xlPos+1) + crossLinkerName + outerPep.subSequence(xlPos+1,outerPep.length());                

                int innerStart = i;

                if (!peptideCanSelfCrosslink) {
                    innerStart = i + 1;
                }

                for (int j = innerStart; j < peptides.length; j++) {
                    String innerPepPlusModString = (String) peptides[j];
                    temp = innerPepPlusModString.split("_");
                    String innerPep = temp[0];
                    String innerPepModString = temp[1];
                    Double innerPepMass = pepMasses.get(innerPepPlusModString);

                    ArrayList<Integer> xlInnerPositions = new ArrayList();

                    pos = innerPep.indexOf(crossLinkedResidue);
                    while (pos != -1) {
                        if ((noCrossLinksOnCTerminus && pos != innerPep.length() - 1) && (noCrossLinksOnModifiedRes && innerPepModString.subSequence(pos + 1, pos + 2).equals("0"))) {
                            xlInnerPositions.add(pos + 1);  //Add position plus one, since the modstring has the n-terminus character at the beginning
                        }
                        pos = innerPep.indexOf(crossLinkedResidue, pos + 1);
                    }

                    for (int y = 0; y < xlInnerPositions.size(); y++) {
                        int xlInnerPos = xlInnerPositions.get(y);
                        String newInnerPepModString = innerPepModString.subSequence(0, xlInnerPos) + "X" + innerPepModString.subSequence(xlInnerPos + 1, innerPepModString.length());
                        String newInnerPepPlusModString = innerPep + "_" + newInnerPepModString;
                        //System.out.println("new inner X:" + newInnerPepPlusModString);                       
                        String crosslinkedPep = outerPepPlusModString + "----" + newInnerPepPlusModString;
                        Double crosslinkedMass = outerPepMass + innerPepMass + crossLinkerMass;
                        pepMasses.put(crosslinkedPep, crosslinkedMass);
                    }
                }
            }
        }


        return pepMasses;

    }

    /*
     * Returns a hashmap of all fragmented ions with their neutral masses
     * HashMap String = peptide sequence (plus special characters for crosslinker)
     * HashMap HashMap String = fragment ion sequence (plus special characters for crosslinker and ion type, e.g. b1, y1.)
     * HashMap HashMap Double = neutral mass of fragment ion
     *
     */
    public HashMap<String, LinkedHashMap<String, Double>> getAllFragmentIonMasses() {

        HashMap<String, LinkedHashMap<String, Double>> allFragmentIonMasses = new HashMap();

        //Todo: it might be better to make allPeptideMasses a global vairable in the class to avoid duplicated function calls
        HashMap<String, Double> allPeptideMasses = getAllPeptideMasses();

        for (String pepPlusModString : allPeptideMasses.keySet()) {

            if (pepPlusModString.contains("----")) {
                //Todo: skip the crosslinking peptide for now
                String[] pepPlusModStrings = pepPlusModString.split("----");

                String pepPlusModString0 = pepPlusModStrings[0];
                Double pepMr0 = getPepMr(pepPlusModString0);


                String pepPlusModString1 = pepPlusModStrings[1];
                Double pepMr1 = getPepMr(pepPlusModString1);

                //getting fragment ions by splitting the frist peptide
                LinkedHashMap<String, Double> tempFragIonMasses = getFragmentIonMasses(pepPlusModString0);
                LinkedHashMap<String, Double> fragIonMasses = new LinkedHashMap();
                for (String pepPlusModPlusIonString : tempFragIonMasses.keySet()) {
                    String[] temp = pepPlusModPlusIonString.split("_");
                    String peptide = temp[0];
                    String modString = temp[1];
                    String ionString = temp[2];
                    Double fragIonMass = tempFragIonMasses.get(pepPlusModPlusIonString);
                    if (modString.contains("X")) {
                        //create a new pepPlusModPlusIonString
                        String newPepPlusModPlusIonString = peptide + "_" + modString + "----" + pepPlusModString1 + "_" + ionString;
                        //get the new fragment ion mass by adding the remain part of the peptide and the other peptide mass and mass of crosslinker
                        Double newFragIonMass = fragIonMass.doubleValue() + pepMr1.doubleValue() + crossLinkerMass;
                        //store in the linkedhashmap
                        fragIonMasses.put(newPepPlusModPlusIonString, newFragIonMass);
                    }
                    else {
                        fragIonMasses.put(pepPlusModPlusIonString, fragIonMass);
                    }
                }

                //getting fragment ions by splitting the second peptide
                tempFragIonMasses = getFragmentIonMasses(pepPlusModString1);
                //LinkedHashMap<String, Double> fragIonMasses1 = new LinkedHashMap();
                for (String pepPlusModPlusIonString : tempFragIonMasses.keySet()) {
                    String[] temp = pepPlusModPlusIonString.split("_");
                    String peptide = temp[0];
                    String modString = temp[1];
                    String ionString = temp[2];
                    Double fragIonMass = tempFragIonMasses.get(pepPlusModPlusIonString);
                    if (modString.contains("X")) {
                        //create a new pepPlusModPlusIonString
                        String newPepPlusModPlusIonString = pepPlusModString0 + "----" + peptide + "_" + modString + "_" + ionString;
                        //get the new fragment ion mass by adding the remain part of the peptide and the other peptide mass and mass of crosslinker
                        Double newFragIonMass = fragIonMass.doubleValue() + pepMr0.doubleValue() + crossLinkerMass;
                        //store in the linkedhashmap
                        fragIonMasses.put(newPepPlusModPlusIonString, newFragIonMass);
                    }
                    else {
                        fragIonMasses.put(pepPlusModPlusIonString, fragIonMass);
                    }
                }


                allFragmentIonMasses.put(pepPlusModString, fragIonMasses);
                //allFragmentIonMasses.put(pepPlusModString, fragIonMasses1);
            }
            else {
                LinkedHashMap<String, Double> fragIonMasses = getFragmentIonMasses(pepPlusModString);
                allFragmentIonMasses.put(pepPlusModString, fragIonMasses);
            }
        }

        return allFragmentIonMasses;
    }

    //Powerset code from http://stackoverflow.com/questions/1670862/obtaining-powerset-of-a-set-in-java
//    public static <T> Set<Set<T>> powerSet(Set<T> originalSet) {
//        Set<Set<T>> sets = new HashSet<Set<T>>();
//        if (originalSet.isEmpty()) {
//            sets.add(new HashSet<T>());
//            return sets;
//        }
//        List<T> list = new ArrayList<T>(originalSet);
//        T head = list.get(0);
//        Set<T> rest = new HashSet<T>(list.subList(1, list.size()));
//        for (Set<T> set : powerSet(rest)) {
//            Set<T> newSet = new HashSet<T>();
//            newSet.add(head);
//            newSet.addAll(set);
//            sets.add(newSet);
//            sets.add(set);
//        }
//        return sets;
//    }

    /*
     * Helper method to create an empty mod string for a given peptide i.e. for PEPTIDER return 0000000000 where each 0 corresponds to unmodified n-terminus (pos 0), unmodified each amino acid and unmodified N-terminuis
     */
    private String makeUnmodifiedModString(String peptide) {

        String modString = "";
        for (int i = 0; i < peptide.length() + 2; i++) {
            modString += 0;
        }
        return modString;
    }

    /*
     * Helper method to insert a mod into a modstring i.e. PMER_000000 to PMER_001000 for Oxidation on M, where oxidation on M has ID = 1
     */
    private String insertModIntoModString(String modString, String varMod,
                                          Set<Integer> modSet) {

        String newModString = modString;

        for (int pos : modSet) {
            pos++;       //need to increment the position, since the modString has N-terminal mod at the beginning
            int modID = varModStringToUniqueID.get(varMod);
            newModString = newModString.subSequence(0, pos) + "" + modID + newModString.subSequence(pos + 1, modString.length());
        }

        return newModString;
    }

    /*
     * Helper method to get all possible fragmented ions (b ion or y ion) from ONE peptides and their neutral masses
     * @param peptidePlusModString the peptide sequence with modification string without crosslink string
     */
    private LinkedHashMap<String, Double> getFragmentIonMasses(
            String peptidePlusModString) {

        LinkedHashMap<String, Double> fragIonMasses = new LinkedHashMap<String, Double>();

        peptidePlusModString = peptidePlusModString.toUpperCase();
        peptidePlusModString = peptidePlusModString.trim();
        //System.out.println("Peptide:" + peptidePlusModString);

        String[] temp = peptidePlusModString.split("_");
        String peptide = temp[0];

        // remove the mod string for C-terminus and N-terminus
        String modString = temp[1].substring(1, temp[1].length() - 1);

        for (int i = 0; i < peptide.length() - 1; i++) {

            // b ion
            String bIon = peptide.substring(0, i + 1);
            String bIonModString = modString.substring(0, i + 1);
            String bIonPlusModString = bIon + "_" + bIonModString;
            Double bIonMass = getFragIonNeutralMass(bIonPlusModString);
            String bIonKey = bIonPlusModString + "_" + "b" + bIon.length();

//            System.out.println("bIonKey=" + bIonKey);

            fragIonMasses.put(bIonKey, bIonMass);

            // y ion
            String yIon = peptide.substring(i + 1, peptide.length());
            String yIonModString = modString.substring(i + 1, modString.length());
            String yIonPlusModString = yIon + "_" + yIonModString;
            Double yIonMass = getFragIonNeutralMass(yIonPlusModString) + CTERM + HMASS;
            String yIonKey = yIonPlusModString + "_" + "y" + yIon.length();

//            System.out.println("yIonKey=" + yIonKey);

            fragIonMasses.put(yIonKey, yIonMass);

        }

        return fragIonMasses;
    }

    /*
     * Helper method to get RESIDUE mass ([M]) of fragment ion with modification
     *
     * The table below is from http://www.matrixscience.com/help/fragmentation_help.html
     *
     * -----------------------------------------------
     * Ion Type | Neutral Mr
     * -------------|---------------------------------
     * a | [N]+[M]-CHO
     * a* | a-NH3
     * a0 | a-H2O
     * b | [N]+[M]-H
     * b* | b-NH3
     * b0 | b-H2O
     * c | [N]+[M]+NH2
     * d | a - partial side chain
     * v | y - complete side chain
     * w | z - partial side chain
     * x | [C]+[M]+CO-H
     * y | [C]+[M]+H
     * y* | y-NH3
     * y0 | y-H2O
     * z | [C]+[M]-NH2
     * ------------------------------------------------
     *
     * [M] is molecular mass of the neutral amino acid residues.
     *
     * [N] is the molecular mass of the neutral N-terminal group
     * [C] is the molecular mass of the neutral C-terminal group
     *
     * To obtain m/z values, add or substract protons as requried to
     * obtain the required charge and divide by the number of charges.
     */
    private double getFragIonNeutralMass(String ionPlusModString) {
        double mass = 0.0;

        String[] temp = ionPlusModString.split("_");
        String ion = temp[0];
        String modString = temp[1];

        // calculate unmodified ion mass 
        String[] aas = ion.split("");
        for (String aa : aas) {
            if (aa.length() == 1) {
                if (aaMap.get(aa) != null) {
                    mass += aaMap.get(aa);
                }
                else {
                    System.out.print("Non-recognized character in fragmented ion:" + ion);
                }
            }
        }

        String[] mods;
        //Code still functional if no modstring is provided
        if (modString != null) {
            mods = modString.split("");
            //Now loop through the var mod string
            for (String mod : mods) {
                if (!mod.equals("0") && !mod.equals("") && !mod.equals("X")) {
                    String varMod = varModUniqueIDToString.get(Integer.parseInt(mod));
                    //System.out.print("mod" + mod + " varMmod:" + varMod);
                    Double varModMass = varMods.get(varMod);
                    //System.out.println(" modmass:" + varModMass);
                    mass += varModMass;
                }
            }
        }

        return mass;
    }

}
